# TODO: Add comment
# 
# Author: broken
###############################################################################
#import library
library(base)
library(raster)
library(gstat)


#path

createKrigingImage = function(regressPm_file){
	
	#PM values
	pmRaster=raster(regressPm_file)
	pm=values(pmRaster)
	corxy=coordinates(pmRaster)
	x=corxy[,'x']
	y=corxy[,'y']
	

	totalCell=length(pmRaster)
	cell = c(1:totalCell)

	table=data.frame(cell,x,y,pm)
	newTable=table
	testTable=subset(table,pm<0)
	trainTable=subset(table,pm>=0)
	
	#caculate variogram
	empiVario=variogram(pm~1,locations=~x+y,data=trainTable)
	
	#sph fit
	sill=min(empiVario$gamma)
	#sphModel=vgm(psill=sill,model="Sph",nugget=0,range=min(empiVario$dist))
	sphModel=vgm(model="Sph",nugget=0,range=1)		
	sphFit=fit.variogram(empiVario,sphModel)
	
	universal_result=krige(id="pm",formula=pm~x+y,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
	
	#edit tiff
	universalPMRaster=pmRaster
	universalPMValue=universal_result[,3]
	universalPMRaster[1:totalCell]=universalPMValue
	
	#edit error tiff
	errorPMRaster=pmRaster
	errorPMValue=universal_result[,4]
	errorPMRaster[1:totalCell]=errorPMValue
	
	#save uk result to tiff
	uk_file = str_replace(regressPm_file,"rg.tif","uk.tif")
	writeRaster(universalPMRaster,filename=uk_file,format="GTiff")
	
	#save uk error to tiff
	error_file = str_replace(regressPm_file,"rg.tif","error.tif")
	writeRaster(errorPMRaster,filename=error_file,format="GTiff")
	print(uk_file)
	
}
createKrigingImage("C:/MOD04_L2.A2014024.0310.051.2014028020029_rg.tif")






